"""
DaoMath Operators (法则数学运算)
基于 ±0 先天均衡的运算体系
"""

import numpy as np
from typing import Union, Callable, List, Tuple
from dataclasses import dataclass
from core.axioms import InnateBalance, DualPair

Number = Union[int, float, complex, np.ndarray]


# ============================================================================
# 核心数据结构
# ============================================================================

@dataclass
class CoherentState:
    """相干态：阴阳相干动态均衡"""
    yin: complex  # 阴分量 -0
    yang: complex  # 阳分量 +0
    coherence: float  # 相干度 [0, 1]
    
    def is_balanced(self, threshold: float = 1e-6) -> bool:
        """检查是否处于均衡态"""
        return abs(self.yin + self.yang) < threshold
    
    def decoherence_risk(self) -> float:
        """退相干风险评估"""
        imbalance = abs(self.yin + self.yang)
        return 1 - np.exp(-imbalance)


@dataclass
class OrderLevel:
    """阶序层级"""
    level: int  # 阶序等级 n
    dimension: int  # 维度 (±0, ±1, ±2, ±3, ...)
    energy: float  # 能量尺度
    coherence_threshold: float  # 相干阈值
    
    def can_transition_to(self, target: 'OrderLevel') -> bool:
        """检查能否跃迁到目标阶序（不可跳跃）"""
        return abs(self.level - target.level) == 1


class DaoOperators:
    """法则数学运算系统"""
    
    # ========================================================================
    # 1. 平衡算子 (Balance Operators) - 基于绝对均衡 -0
    # ========================================================================
    
    @staticmethod
    def balance_sum(a: Number, b: Number) -> Tuple[Number, Number, Number]:
        """
        均衡加法: a ⊕ b → {正分量, 负分量, 中性分量}
        
        数学原理：
        任意数可分解为 ±0 基底：x = x⁺ + x⁻ + x⁰
        其中 x⁺ ≥ 0, x⁻ ≤ 0, x⁺ + x⁻ + x⁰ = x
        """
        total = a + b
        
        if isinstance(total, complex):
            positive = max(0, total.real) + max(0, total.imag) * 1j
            negative = min(0, total.real) + min(0, total.imag) * 1j
            neutral = 0 + 0j
        else:
            positive = max(0, float(total))
            negative = min(0, float(total))
            neutral = 0.0
        
        # 验证守恒
        assert abs(positive + negative + neutral - total) < 1e-10
        
        return positive, negative, neutral
    
    @staticmethod
    def balance_product(a: Number, b: Number) -> InnateBalance:
        """
        均衡乘法: a ⊗ b → ±0
        
        符号规则（阴阳法则）：
        + × + = + (阳生阳)
        - × - = + (阴生阳)
        + × - = - (阴阳相克)
        """
        product = a * b
        
        if isinstance(product, complex):
            # 复数分解为实部（阳）和虚部（阴）
            real_sign = np.sign(product.real) if product.real != 0 else 0
            imag_sign = np.sign(product.imag) if product.imag != 0 else 0
            
            return InnateBalance(
                positive=complex(max(0, product.real), max(0, product.imag)),
                negative=complex(min(0, product.real), min(0, product.imag)),
                innate_zero=0j
            )
        else:
            if product > 0:
                return InnateBalance(positive=product, negative=0, innate_zero=0)
            elif product < 0:
                return InnateBalance(positive=0, negative=product, innate_zero=0)
            else:
                return InnateBalance(positive=0, negative=0, innate_zero=0)
    
    @staticmethod
    def annihilation(yin: Number, yang: Number) -> Number:
        """
        湮灭归零: yin ⊖ yang → 0
        
        物理意义：正负粒子对湮灭，能量守恒转化
        数学形式：-x + (+x) = 0
        """
        return yin + yang
    
    # ========================================================================
    # 2. 对偶算子 (Duality Operators) - 互为镜像对称
    # ========================================================================
    
    @staticmethod
    def dual(x: Number) -> Number:
        """
        对偶算子: † 
        镜像对称：†x = -x
        """
        return -x
    
    @staticmethod
    def conjugate_dual(x: complex) -> complex:
        """
        共轭对偶: * 
        复数镜像：z* = a - bi (若 z = a + bi)
        """
        return np.conj(x)
    
    @staticmethod
    def inverse_dual(x: Number) -> Number:
        """
        逆对偶: ⁻¹
        乘法镜像：x⁻¹ · x = 1
        """
        if x == 0:
            raise ValueError("Cannot invert zero (违反守恒)")
        return 1 / x
    
    @staticmethod
    def create_dual_pair(x: Number) -> DualPair[Number]:
        """创建对偶对：{x, †x}"""
        return DualPair(positive=x, negative=-x)
    
    # ========================================================================
    # 3. 三元算子 (Ternary Operators) - 生克制化
    # ========================================================================
    
    @staticmethod
    def ternary_mean(a: Number, b: Number, c: Number) -> Number:
        """三元均值：均衡态"""
        return (a + b + c) / 3
    
    @staticmethod
    def ternary_product(a: Number, b: Number, c: Number) -> Number:
        """三元乘积：守恒量"""
        return a * b * c
    
    @staticmethod
    def ternary_balance(positive: Number, negative: Number, neutral: Number) -> Number:
        """
        三元平衡态合成
        ±0 = {+1, -1, 0} → 总和为0
        """
        result = positive + negative + neutral
        return result
    
    @staticmethod
    def wuxing_cycle(elements: List[Number]) -> List[Number]:
        """
        五行相生相克循环
        
        相生：木→火→土→金→水→木
        相克：木→土→水→火→金→木
        """
        if len(elements) != 5:
            raise ValueError("五行需要5个元素")
        
        # 相生（顺时针）
        sheng = [
            (elements[0] + elements[1]) / 2,  # 木生火
            (elements[1] + elements[2]) / 2,  # 火生土
            (elements[2] + elements[3]) / 2,  # 土生金
            (elements[3] + elements[4]) / 2,  # 金生水
            (elements[4] + elements[0]) / 2,  # 水生木
        ]
        
        return sheng
    
    # ========================================================================
    # 4. 守恒算子 (Conservation Operators) - 绝对守恒 +0
    # ========================================================================
    
    @staticmethod
    def conserved_sum(values: List[Number]) -> Number:
        """
        守恒和: Σᵢ vᵢ = const
        
        物理意义：能量守恒、电荷守恒
        """
        return sum(values)
    
    @staticmethod
    def conserved_product(values: List[Number]) -> Number:
        """
        守恒积: Πᵢ vᵢ = const
        
        物理意义：角动量守恒
        """
        result = 1
        for v in values:
            result *= v
        return result
    
    @staticmethod
    def phase_space_volume(positions: np.ndarray, momenta: np.ndarray) -> float:
        """
        相空间体积 (Liouville定理)
        
        物理意义：Hamilton系统中，相空间体积守恒
        数学形式：∫∫ dq dp = const
        """
        if len(positions) != len(momenta):
            raise ValueError("位置和动量维度必须相同")
        
        # 计算 Π(Δqᵢ · Δpᵢ)
        volume = 1.0
        for q, p in zip(positions, momenta):
            volume *= abs(q * p)
        
        return float(volume)
    
    @staticmethod
    def noether_charge(symmetry_generator: Callable, states: List[Number]) -> Number:
        """
        Noether守恒荷
        
        对称性 → 守恒律（Noether定理）
        时间平移对称 → 能量守恒
        空间平移对称 → 动量守恒
        旋转对称 → 角动量守恒
        """
        charges = [symmetry_generator(state) for state in states]
        return sum(charges)
    
    # ========================================================================
    # 5. 涌现算子 (Emergence Operators) - 分形递归
    # ========================================================================
    
    @staticmethod
    def emergent_mean(lower_level: List[Number]) -> Number:
        """从下层涌现均值（一阶统计量）"""
        return np.mean(lower_level)
    
    @staticmethod
    def emergent_variance(lower_level: List[Number]) -> float:
        """从下层涌现方差（二阶统计量）"""
        return float(np.var(lower_level))
    
    @staticmethod
    def emergent_entropy(probabilities: List[float]) -> float:
        """
        从微观态涌现熵（信息论）
        
        Shannon熵：S = -Σᵢ pᵢ log pᵢ
        物理意义：系统的不确定性/信息量
        """
        probs = np.array(probabilities)
        probs = probs[probs > 0]  # 避免log(0)
        return float(-np.sum(probs * np.log2(probs)))
    
    @staticmethod
    def fractal_dimension(data: np.ndarray, epsilon: float = 0.01) -> float:
        """
        分形维数（自相似性度量）
        
        数学形式：D = lim(ε→0) log N(ε) / log(1/ε)
        应用：海岸线长度、股票波动
        """
        scales = np.logspace(-2, 0, 20)
        counts = []
        
        for scale in scales:
            boxes = len(np.unique((data / scale).astype(int)))
            counts.append(boxes)
        
        # 线性拟合 log-log 图
        coeffs = np.polyfit(np.log(scales), np.log(counts), 1)
        return float(-coeffs[0])  # 斜率的负值
    
    @staticmethod
    def recursive_emergence(
        base: Number,
        rule: Callable[[Number], Number],
        levels: int
    ) -> List[Number]:
        """
        递归涌现：从低阶序到高阶序
        
        n=0: 基础态
        n=1: rule(n=0)
        n=2: rule(rule(n=0))
        ...
        """
        results = [base]
        current = base
        
        for _ in range(levels):
            current = rule(current)
            results.append(current)
        
        return results
    
    # ========================================================================
    # 6. 阶序算子 (Order Operators) - 阶序不可跳跃
    # ========================================================================
    
    @staticmethod
    def order_transition(
        current: OrderLevel,
        target: OrderLevel,
        energy_input: float
    ) -> Tuple[bool, str]:
        """
        阶序跃迁检验
        
        规则：
        1. 只能相邻跃迁（n → n±1）
        2. 需要足够能量（熵减需要能量输入）
        3. 向上跃迁：熵减（需要外部能量）
        4. 向下跃迁：熵增（自发过程）
        """
        if not current.can_transition_to(target):
            return False, "违反阶序不可跳跃原则"
        
        delta_level = target.level - current.level
        
        if delta_level > 0:  # 向上跃迁（熵减）
            required_energy = target.energy - current.energy
            if energy_input < required_energy:
                return False, f"能量不足：需要 {required_energy}，提供 {energy_input}"
            return True, "熵减跃迁成功"
        
        else:  # 向下跃迁（熵增）
            return True, "退相干衰减（自发过程）"
    
    @staticmethod
    def entropy_production(
        initial_order: OrderLevel,
        final_order: OrderLevel
    ) -> float:
        """
        熵产生计算
        
        ΔS = S_final - S_initial
        ΔS > 0: 熵增（退相干）
        ΔS < 0: 熵减（需要外部能量）
        ΔS = 0: 可逆过程
        """
        # 简化模型：熵与阶序能量成反比
        S_initial = 1.0 / (initial_order.energy + 1e-10)
        S_final = 1.0 / (final_order.energy + 1e-10)
        
        return S_final - S_initial
    
    # ========================================================================
    # 7. 相干算子 (Coherence Operators) - 阴阳相干/退相干
    # ========================================================================
    
    @staticmethod
    def coherence_measure(yin: complex, yang: complex) -> float:
        """
        相干度量
        
        C = |⟨yin|yang⟩| / (||yin|| · ||yang||)
        C = 1: 完全相干
        C = 0: 完全退相干
        """
        if abs(yin) < 1e-10 or abs(yang) < 1e-10:
            return 0.0
        
        inner_product = np.conj(yin) * yang
        norm_product = abs(yin) * abs(yang)
        
        return float(abs(inner_product) / norm_product)
    
    @staticmethod
    def decoherence_rate(
        state: CoherentState,
        environment_coupling: float,
        temperature: float
    ) -> float:
        """
        退相干率
        
        Γ = γ · (1 + n̄)
        γ: 环境耦合强度
        n̄: 热激发数（温度相关）
        
        物理意义：系统与环境相互作用导致相干性丧失
        """
        # Bose-Einstein分布
        k_B = 1.380649e-23  # Boltzmann常数
        hbar = 1.054571817e-34  # 约化Planck常数
        omega = 1e15  # 特征频率（Hz）
        
        if temperature > 0:
            n_bar = 1.0 / (np.exp(hbar * omega / (k_B * temperature)) - 1)
        else:
            n_bar = 0
        
        gamma = environment_coupling * (1 + n_bar)
        return float(gamma)
    
    @staticmethod
    def coherence_revival(
        decohered_state: CoherentState,
        control_energy: float
    ) -> CoherentState:
        """
        相干恢复（量子纠错）
        
        通过外部控制恢复部分相干性
        """
        revival_factor = min(1.0, control_energy / 10.0)  # 简化模型
        
        new_coherence = decohered_state.coherence + revival_factor * (1 - decohered_state.coherence)
        
        return CoherentState(
            yin=decohered_state.yin,
            yang=decohered_state.yang,
            coherence=new_coherence
        )
    
    # ========================================================================
    # 8. 循环算子 (Cyclic Operators) - 因果闭环
    # ========================================================================
    
    @staticmethod
    def causal_loop(
        initial_state: Number,
        forward_map: Callable[[Number], Number],
        backward_map: Callable[[Number], Number],
        iterations: int
    ) -> Tuple[List[Number], bool]:
        """
        因果闭环检验
        
        检查 backward(forward(x)) ≈ x（自洽性）
        """
        states = [initial_state]
        current = initial_state
        
        for _ in range(iterations):
            # 正向演化
            forward_state = forward_map(current)
            # 反向演化
            backward_state = backward_map(forward_state)
            
            states.append(backward_state)
            current = backward_state
        
        # 检查是否回到初态
        is_closed = abs(states[-1] - initial_state) < 1e-6
        
        return states, is_closed
    
    @staticmethod
    def spiral_evolution(
        center: complex,
        radius: float,
        angular_velocity: float,
        time: np.ndarray
    ) -> np.ndarray:
        """
        螺旋演化（0既是起点也是归0点）
        
        z(t) = center + r(t) · exp(iωt)
        r(t) = r₀ · exp(-αt)  (阻尼螺旋)
        """
        damping = 0.1
        r_t = radius * np.exp(-damping * time)
        phase = angular_velocity * time
        
        trajectory = center + r_t * np.exp(1j * phase)
        
        return trajectory
    
    # ========================================================================
    # 9. 阈值算子 (Threshold Operators) - 均衡上下限
    # ========================================================================
    
    @staticmethod
    def threshold_check(
        value: Number,
        lower_bound: Number,
        upper_bound: Number
    ) -> Tuple[bool, str]:
        """
        阈值检查（红线机制）
        
        失衡超出阈值 → 触发守恒归零
        """
        if value < lower_bound:
            return False, f"突破下限阈值 {lower_bound}，触发归零"
        elif value > upper_bound:
            return False, f"突破上限阈值 {upper_bound}，触发归零"
        else:
            return True, "在均衡阈值内"
    
    @staticmethod
    def adaptive_threshold(
        history: List[Number],
        sensitivity: float = 3.0
    ) -> Tuple[Number, Number]:
        """
        自适应阈值
        
        基于历史数据动态调整阈值
        μ ± sensitivity · σ
        """
        mean = np.mean(history)
        std = np.std(history)
        
        lower = mean - sensitivity * std
        upper = mean + sensitivity * std
        
        return lower, upper
    
    # ========================================================================
    # 10. 维度算子 (Dimension Operators) - 从低维到高维
    # ========================================================================
    
    @staticmethod
    def dimension_embedding(
        low_dim_data: np.ndarray,
        target_dim: int,
        method: str = 'zero_padding'
    ) -> np.ndarray:
        """
        维度嵌入：从低维到高维
        
        ±0维 → ±1维 → ±2维 → ±3维
        """
        current_dim = low_dim_data.shape[-1] if low_dim_data.ndim > 1 else 1
        
        if target_dim < current_dim:
            raise ValueError("不能降维（违反阶序原则）")
        
        if method == 'zero_padding':
            padding = np.zeros(target_dim - current_dim)
            return np.concatenate([low_dim_data.flatten(), padding])
        
        elif method == 'fractal':
            # 分形扩展
            result = low_dim_data.flatten()
            while len(result) < target_dim:
                result = np.concatenate([result, result * 0.5])
            return result[:target_dim]
        
        else:
            raise ValueError(f"未知方法: {method}")
    
    @staticmethod
    def dimension_projection(
        high_dim_data: np.ndarray,
        target_dim: int,
        method: str = 'pca'
    ) -> np.ndarray:
        """
        维度投影：从高维到低维（退相干）
        
        注意：这是有损过程（熵增）
        """
        if method == 'pca':
            # 简化PCA
            from sklearn.decomposition import PCA
            pca = PCA(n_components=target_dim)
            return pca.fit_transform(high_dim_data.reshape(1, -1))[0]
        
        elif method == 'truncation':
            return high_dim_data[:target_dim]
        
        else:
            raise ValueError(f"未知方法: {method}")


# ============================================================================
# 完整示例
# ============================================================================

if __name__ == "__main__":
    ops = DaoOperators()
    
    print("=" * 70)
    print("法则数学运算系统完整示例")
    print("=" * 70)
    print()
    
    # ========== 1. 平衡算子 ==========
    print("【1. 平衡算子 - 绝对均衡 -0】")
    print("-" * 70)
    a, b = 3, -5
    pos, neg, neu = ops.balance_sum(a, b)
    print(f"均衡加法: {a} ⊕ {b}")
    print(f"  正分量: {pos}")
    print(f"  负分量: {neg}")
    print(f"  中性分量: {neu}")
    print(f"  验证守恒: {pos + neg + neu} = {a + b}")
    
    # 湮灭归零
    yin, yang = -5, 5
    result = ops.annihilation(yin, yang)
    print(f"\n湮灭归零: {yin} ⊖ {yang} = {result}")
    print()
    
    # ========== 2. 对偶算子 ==========
    print("【2. 对偶算子 - 镜像对称】")
    print("-" * 70)
    x = 3 + 4j
    print(f"原值: x = {x}")
    print(f"对偶: †x = {ops.dual(x)}")
    print(f"共轭: x* = {ops.conjugate_dual(x)}")
    print(f"逆对偶: x⁻¹ = {ops.inverse_dual(x):.4f}")
    
    pair = ops.create_dual_pair(5)
    print(f"\n对偶对: {pair}")
    print()
    
    # ========== 3. 三元/五行算子 ==========
    print("【3. 三元算子 - 生克制化】")
    print("-" * 70)
    elements = [1.0, 2.0, 3.0, 4.0, 5.0]  # 木火土金水
    cycle = ops.wuxing_cycle(elements)
    print("五行相生:")
    labels = ["木生火", "火生土", "土生金", "金生水", "水生木"]
    for label, value in zip(labels, cycle):
        print(f"  {label}: {value:.2f}")
    print()
    
    # ========== 4. 守恒算子 ==========
    print("【4. 守恒算子 - 绝对守恒 +0】")
    print("-" * 70)
    values = [1, 2, 3, 4, 5]
    print(f"守恒和: Σ{values} = {ops.conserved_sum(values)}")
    print(f"守恒积: Π{values} = {ops.conserved_product(values)}")
    
    # 相空间体积
    q = np.array([1.0, 2.0, 3.0])
    p = np.array([0.5, 1.0, 1.5])
    volume = ops.phase_space_volume(q, p)
    print(f"相空间体积: {volume:.4f}")
    print()
    
    # ========== 5. 涌现算子 ==========
    print("【5. 涌现算子 - 分形递归】")
    print("-" * 70)
    data = [1, 2, 3, 4, 5]
    print(f"涌现均值: {ops.emergent_mean(data):.2f}")
    print(f"涌现方差: {ops.emergent_variance(data):.2f}")
    
    probs = [0.2, 0.3, 0.5]
    entropy = ops.emergent_entropy(probs)
    print(f"Shannon熵: {entropy:.4f} bits")
    
    # 递归涌现
    recursive = ops.recursive_emergence(1, lambda x: x * 2, 5)
    print(f"递归涌现: {recursive}")
    print()
    
    # ========== 6. 阶序算子 ==========
    print("【6. 阶序算子 - 不可跳跃】")
    print("-" * 70)
    level_0 = OrderLevel(0, 0, 1.0, 0.9)  # 基态
    level_1 = OrderLevel(1, 1, 10.0, 0.8)  # 一阶激发态
    level_2 = OrderLevel(2, 2, 50.0, 0.7)  # 二阶激发态
    
    # 尝试跃迁
    success, msg = ops.order_transition(level_0, level_1, 12.0)
    print(f"0→1 跃迁 (能量12): {msg}")
    
    success, msg = ops.order_transition(level_0, level_2, 100.0)
    print(f"0→2 跃迁 (能量100): {msg}")
    
    # 熵产生
    delta_S = ops.entropy_production(level_0, level_1)
    print(f"熵变化: ΔS = {delta_S:.4f}")
    print()
    
    # ========== 7. 相干算子 ==========
    print("【7. 相干算子 - 阴阳相干】")
    print("-" * 70)
    yin = 1 + 0j
    yang = 0 + 1j
    coherence = ops.coherence_measure(yin, yang)
    print(f"相干度: C = {coherence:.4f}")
    
    state = CoherentState(yin, yang, 0.9)
    gamma = ops.decoherence_rate(state, 0.1, 300)
    print(f"退相干率: Γ = {gamma:.2e} Hz")
    print()
    
    # ========== 8. 循环算子 ==========
    print("【8. 循环算子 - 因果闭环】")
    print("-" * 70)
    states, closed = ops.causal_loop(
        1.0,
        lambda x: x * 2,
        lambda x: x / 2,
        3
    )
    print(f"因果闭环: {states}")
    print(f"自洽性: {'通过' if closed else '失败'}")
    
    # 螺旋演化
    t = np.linspace(0, 10, 100)
    trajectory = ops.spiral_evolution(0+0j, 1.0, 1.0, t)
    print(f"螺旋演化: 起点 {trajectory[0]:.4f} → 终点 {trajectory[-1]:.4f}")
    print()
    
    # ========== 9. 阈值算子 ==========
    print("【9. 阈值算子 - 红线机制】")
    print("-" * 70)
    history = [1, 2, 3, 4, 5]
    lower, upper = ops.adaptive_threshold(history, sensitivity=2.0)
    print(f"自适应阈值: [{lower:.2f}, {upper:.2f}]")
    
    test_values = [3, 0, 10]
    for val in test_values:
        ok, msg = ops.threshold_check(val, lower, upper)
        status = "✓" if ok else "✗"
        print(f"  {status} 测试值 {val}: {msg}")
    print()
    
    # ========== 10. 维度算子 ==========
    print("【10. 维度算子 - 阶序演化】")
    print("-" * 70)
    data_1d = np.array([1, 2, 3])
    data_3d = ops.dimension_embedding(data_1d, 6, method='fractal')
    print(f"1维 → 6维嵌入: {data_3d}")
    print()
    
    print("=" * 70)
    print("示例完成")
    print("=" * 70)
